#### Prepares results from model-single.R for model-combine.R.

## NOTE: Assumes that output.suffix, allfits, sumstats already defined
## Creates: mcmc, regions, uniregs

names(sumstats)[1] <- "fullparam"
sumstats$param <- gsub("\\[\\d+\\]", "", gsub("_coeff\\[(\\d+)\\]", "_coeff\\1", sumstats$fullparam))
names(sumstats)[names(sumstats) == 'n_eff'] <- 'neff' # if used old convention

sumstats <- subset(sumstats, region %in% unique(allfits$region)) # in case working with subset

mcmc <- allfits[, -1]
regions <- allfits$region
uniregs <- unique(regions)

## Drop outlier municipalities
for (cc in 1:ncol(mcmc))
    mcmc[!is.finite(mcmc[, cc]), cc] <- NA
if (length(grep('-eq24-orig', output.suffix)) > 0 || length(grep('-nos-orig', output.suffix)) > 0 || length(grep('-eq24rational-orig', output.suffix)) > 0 || length(grep('-nosrational-orig', output.suffix)) > 0) {
    mcmc$cost_to_planting[abs(mcmc$cost_to_planting) > quantile(abs(mcmc$cost_to_planting), .95)] <- NA # Adjust after have more values
    mcmc$death_coeff1[abs(mcmc$death_coeff1) > quantile(abs(mcmc$death_coeff1), .95, na.rm=T)] <- NA # Adjust after have more values
    mcmc$cost_to_planting[mcmc$cost_to_planting < quantile(mcmc$cost_to_planting, .05, na.rm=T)] <- NA
}
## Drop unconstrained parameters
for (rr in 1:length(uniregs)) {
    ## Check if we had any EDDs
    if (!all(df$edd1000[df$Município == uniregs[rr]] == 0, na.rm=T))
        next

    ## Drop these parameter estimates
    for (cc in which(names(mcmc) %in% c(paste0('yield_coeff', which(weatherpreds == 'edd1000')), 'death_coeff1')))
        mcmc[regions == uniregs[rr], cc] <- NA
}
## Normalize costs
for (rr in 1:length(uniregs)) {
    denom <- mean(df$total.harvest[df$Município == uniregs[rr]], na.rm=T)
    if (!is.na(denom) && denom == 0)
        denom <- NA
    mcmc$cost_to_planting[regions == uniregs[rr]] <- mcmc$cost_to_planting[regions == uniregs[rr]] / denom
}

## Adjust means and sds within sumstats based on these changes
for (region in uniregs) {
    print(region)
    for (cc in 1:ncol(mcmc)) {
        sumstats$mean[sumstats$param == names(mcmc)[cc] & sumstats$region == region] <- mean(mcmc[regions == region, cc])
        sumstats$sd[sumstats$param == names(mcmc)[cc] & sumstats$region == region] <- sd(mcmc[regions == region, cc])
    }
}

